We will use average annual sea surface temperature (SST) from the years 2008 to 2012 to characterize the average sea surface temperature within the region. The data we are working with was originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.
We will use average annual sea surface temperature (SST) from the years 2008 to 2012 to characterize the average sea surface temperature within the region. The data we are working with was originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.
here() starts at C:/Users/kiran/Documents/MEDS 2022-2023/kiranfavre.github.io
Show code
library(terra)library(tmap)library(tmaptools)#file path set by Rprojsetwd(here())
Show code
#read in West coast EEZ shape file w terra wc_EEZ_regions <-st_read(here("posts","2023-05-03_geospatial","data","wc_regions_clean.shp"))
Reading layer `wc_regions_clean' from data source
`C:\Users\kiran\Documents\MEDS 2022-2023\kiranfavre.github.io\posts\2023-05-03_geospatial\data\wc_regions_clean.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
Geodetic CRS: WGS 84
Show code
#read in SST rastersavg_sst_2008 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2008.tif"))avg_sst_2009 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2009.tif"))avg_sst_2010 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2010.tif"))avg_sst_2011 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2011.tif"))avg_sst_2012 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2012.tif"))#combine into raster stack all_sst <-c(avg_sst_2008, avg_sst_2009, avg_sst_2010, avg_sst_2011, avg_sst_2012)#read in bathymetry rasterdepth <-rast(here("posts","2023-05-03_geospatial","data","depth.tif"))
Checking if the coordinate reference systems are the same. The sea surface temperature data needs to be reprojected to match the coordinate reference systems of the EEZ and depth data.
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Show code
#re project using the terra wayall_sst_reproj <-project(all_sst, wc_EEZ_regions)
Warning: [project,SpatRaster] argument y (the crs) should be a character value
Process data
Next, I need process the SST and depth data so that they can be combined, I will resample to match the SST data using the nearest neighbor approach since the extents, resolutions, and positions of the depth and SST are different.
Show code
#find the mean SST from 2008-2012mean_sst <-mean(all_sst_reproj)#convert sst from K to Cmean_sst_c <- mean_sst -273.15
Show code
#crop depth rast to match the extent of the SST rastdepth_cropped <-crop(depth,mean_sst_c)#using nearest neighbor approach to resampledepth_resampled <-resample(depth_cropped, mean_sst_c,method ="near")#stack to check if they have the same resolutionresolution_test <-c(depth_resampled, all_sst_reproj)#they stacked properly, indicating that they both have same resolutions, extent, and crs
Find suitable locations
To determine suitable locations for marine aquaculture, I need to look for locations that are suitable in terms of both SST and depth.
Show code
#oyster happy spots#sea surface temperature: 11-30°C\#depth: 0-70 meters below sea level#reclassify sst into suitable locations for oysters#reclassification matrix for suitable range of sstrcl_sst <-matrix(c(-Inf, 11, NA, 11, 30,1,30, Inf, NA), ncol =3, byrow =TRUE) #anything outside of range is NA#reclassifying raster using a reclassification matrixsst_suitable_locs <-classify(mean_sst_c,rcl = rcl_sst, include.lowest =TRUE)
Show code
#reclassify depth into suitable locations for oysters#reclassification matrix for suitable range of depthsrcl_depth <-matrix(c(-Inf, -70, NA,-70, 0, 1,0, Inf, NA),ncol =3, byrow =TRUE)#reclassifying raster using a reclassification matrixdepth_suitable_locs <-classify(depth_resampled,rcl = rcl_depth,include.lowest =TRUE)
Show code
#define functionmult_fun <-function(x, y) {return(x*y)}#find locations that satisfy both conditionsoyst_suitable_habitat <-lapp(c(sst_suitable_locs, depth_suitable_locs), mult_fun)
Determining the most suitable EEZ
I want to determine the total suitable area within each EEZ so I can rank zones by priority. To do this, I will need to find the total area of suitable locations within each EEZ.
Show code
# making obj of individual cell size of areas that are suitable for oysters so we can find the size of each cellgrid_cell_size <-cellSize(oyst_suitable_habitat,mask =TRUE, #keep values within the suitable habitattransform =TRUE,unit ="km")plot(grid_cell_size, main ="Cell Size")
Show code
# rasterize wc shape file to help find total suitable area# Rasterize(): Transfer values associated with 'object' type spatial data (points, lines, polygons) to raster cells.# wc_EEZ regions are the points, and we will transfer those values to the raster i set up earlier for suitable oyster habitatswc_rasterized <-rasterize(wc_EEZ_regions, oyst_suitable_habitat,field ="rgn")
Show code
#make a mask of wc raster and suitable locs for oysters#goal is to have areas of suitable oyster habitat that also fall in w coast EEZ regionswc_mask <-mask(wc_rasterized, oyst_suitable_habitat,updatevalue =NA,inverse =FALSE)
Show code
#find area of each suitable zones with zonal() in each of 5 regions of EEZwc_suitable_area <-zonal(grid_cell_size, wc_mask, na.rm =TRUE, sum)
Show code
#join datawc_suitable_EEZ <-full_join(wc_EEZ_regions, wc_suitable_area, by ="rgn") |>mutate(suitable_area = area, #rename area to suitable area percentage_suitable = (suitable_area/area_km2 *100), #percentage of suitable area.before = geometry)print(wc_suitable_EEZ)
Simple feature collection with 5 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
Geodetic CRS: WGS 84
rgn rgn_key area_m2 rgn_id area_km2 area
1 Oregon OR 179994061293 1 179994.06 1074.2562
2 Northern California CA-N 164378809215 2 164378.81 178.0246
3 Central California CA-C 202738329147 3 202738.33 4069.5671
4 Southern California CA-S 206860777840 4 206860.78 3508.1870
5 Washington WA 66898309678 5 66898.31 2378.2758
suitable_area percentage_suitable geometry
1 1074.2562 0.5968287 MULTIPOLYGON (((-123.4318 4...
2 178.0246 0.1083014 MULTIPOLYGON (((-124.2102 4...
3 4069.5671 2.0073003 MULTIPOLYGON (((-122.9928 3...
4 3508.1870 1.6959169 MULTIPOLYGON (((-120.6505 3...
5 2378.2758 3.5550611 MULTIPOLYGON (((-122.7675 4...
Visualize results
Show code
#set to interactive modetmap_mode("view")
tmap mode set to interactive viewing
Show code
#map for total suitable area for oysters by regionoyst_area_map <-tm_basemap("Stamen.Terrain") +tm_shape(wc_suitable_EEZ) +tm_polygons(col ='area',palette ='Oranges',alpha =0.75,border.col ='black',title ="Total Suitable Area") +tm_text("rgn", size =0.54) +tm_scale_bar(position =c("left", "right"))oyst_area_map
Show code
#map for percent suitable area by regionoyst_percentage_map <-tm_basemap("Stamen.Terrain") +tm_shape(wc_suitable_EEZ) +tm_polygons(col ='percentage_suitable',palette ='Purples',alpha =0.75,border.col ='black',title ="Percentage of Suitable Area") +tm_text("rgn", size =0.54) +tm_scale_bar(position =c("left", "right"))oyst_percentage_map
Expanding this workflow to other species
Show code
find_suitable_locs <-function(SST_low, SST_high, depth_low, depth_high, spp_name) {#read in West coast EEZ shape file w terra wc_EEZ_regions <-st_read(here("posts","2023-05-03_geospatial","data","wc_regions_clean.shp"))#read in SST rastersavg_sst_2008 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2008.tif"))avg_sst_2009 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2009.tif"))avg_sst_2010 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2010.tif"))avg_sst_2011 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2011.tif"))avg_sst_2012 <-rast(here("posts","2023-05-03_geospatial","data","average_annual_sst_2012.tif"))#combine into raster stack all_sst <-c(avg_sst_2008, avg_sst_2009, avg_sst_2010, avg_sst_2011, avg_sst_2012)#read in bathymetry rasterdepth <-rast(here("posts","2023-05-03_geospatial","data","depth.tif"))#re project using the terra way all_sst_reproj <-project(all_sst, depth)#find the mean SST from 2008-2012 mean_sst <-mean(all_sst_reproj)#convert sst from K to C mean_sst_c <- mean_sst -273.15#crop depth rast to match the extent of the SST rast depth_cropped <-crop(depth,mean_sst_c)#using nearest neighbor approach to resample depth_resampled <-resample(depth_cropped, mean_sst_c,method ="near")#make reclassification matrix rcl_sst <-matrix(c(-Inf, SST_low, NA, SST_low, SST_high,1, SST_high, Inf, NA), ncol =3, byrow =TRUE) #reclassifying sst raster using a reclassification matrix sst_suitable_locs <-classify(mean_sst_c,rcl = rcl_sst, include.lowest =TRUE)#reclassifying depth raster using a reclassification matrix rcl_depth <-matrix(c(-Inf, depth_low, NA, depth_low, depth_high, 1, depth_high, Inf, NA),ncol =3, byrow =TRUE) depth_suitable_locs <-classify(depth_resampled,rcl = rcl_depth,include.lowest =TRUE)#define function to multiply layers mult_fun <-function(x, y) {return(x*y)}#find locations that satisfy both conditions oyst_suitable_habitat <-lapp(c(sst_suitable_locs, depth_suitable_locs), mult_fun)#find grid cell size grid_cell_size <-cellSize(oyst_suitable_habitat,mask =TRUE, transform =TRUE,unit ="km") wc_rasterized <-rasterize(wc_EEZ_regions, oyst_suitable_habitat,field ="rgn") wc_mask <-mask(wc_rasterized, oyst_suitable_habitat,updatevalue =NA,inverse =FALSE) wc_suitable_area <-zonal(grid_cell_size, wc_mask, na.rm =TRUE, sum) wc_suitable_EEZ <-full_join(wc_EEZ_regions, wc_suitable_area, by ="rgn") |>mutate(suitable_area = area, #rename area to suitable area percentage_suitable = (suitable_area/area_km2 *100), #percentage of suitable area.before = geometry)##make maps #set to interactive viewing tmap_mode("view")#total area viz_total_area <-tm_basemap("Stamen.Terrain")+tm_shape(wc_suitable_EEZ) +tm_polygons(col ='area',palette ='Oranges',alpha =0.75,border.col ='black',title =paste0("Total Suitable Area for ", spp_name ," in West Coast EEZs")) +tm_text("rgn", size =0.54) viz_total_area#percentage suitable viz_percentage <-tm_basemap("Stamen.Terrain") +tm_shape(wc_suitable_EEZ) +tm_polygons(col ='percentage_suitable',palette ='Purples',alpha =0.75,border.col ='black',title =paste0("Percentage of Suitable Area for ", spp_name ," in West Coast EEZs")) +tm_text("rgn", size =0.54) viz_percentagetmap_arrange(viz_total_area, viz_percentage)}
I chose the Red Abalone, Haliotis rufescens. https://www.sealifebase.ca/summary/Haliotis-rufescens.html